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Abstract The rarity of limestone cave species due to habitat degradation makes them of special interest in conservation 
biology. The wild Chinese Giant Salamander Andrias davidianus, an evolutionarily distinct and globally endangered 
species, are nearly all obligate into living in inaccessible mountain caves now. We detected that only 14—29 breeders, 
with the effective population size of 9—25 (5—44 of 95% confidence interval), were in each of three caves, through 
genotyping 20 microsatellite loci on larvae that were flushed out of caves. Both heterzygosity excess and M (ratios of 
allele number to allele size range) tests indicated severe genetic bottlenecks among populations. Both mitochondrial, 
with only one or two haplotypes of D-loop region (770—771Ыр) in each population, and nuclear genetic structure 
showed clear divergence between populations. Considering the long life history, small population size, and genetic 
differentiationof the Chinese Giant Salamander, putting an end to poaching and recovering the karst ecosystem instead 


of releasing may be the sole measures to save this severely threated species. 
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1. Introduction 


Subterranean environments are renowned as fragile 
and vulnerable ecosystems. Anthropogenic activities 
such as agriculture, urban development, waste disposal, 
deforestation, and resource extraction constitute 
environmental pressures leading to polluted aquifers, 
destroyed caves, denuded landscapes, and declining 
biodiversity (van Beynen et al., 2012). Karst caves serve 
as shelters for a variety of plants and animals, some 
of which are endemic to caves and are referred to as 
troglobites (Culver and Pipan, 2009). Among vertebrates, 
few amphibians are restricted to karst caves, such as 
cave-dwelling salamanders Proteus anguinus, Eurycea 
spp., and Necturus sp. (Goricki et al., 2012), but more 
typically, they utilize caves as one of a variety of suitable 
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habitats for reproduction, refuge, and/or hibernation 
(Sket, 1997; Liu et al., 2009; Niemiller and Miller, 2009; 
Kohler et al., 2010; Liu, 2010). The rarity of cave species 
makes them of special interest in conservation biology. 
The Chinese Giant Salamander (CGS) Andrias 
davidianus is the largest extant amphibian with a 
generation length of about 15 years. Historically, this 
species occurred widely in the watersheds of the Yellow, 
Yangtze, and Pearl River. However, it has experienced 
dramatic population declines in the past sixty years, 
estimated to be as much as 80% over the last three 
generations due to over-exploitation for the commercial 
luxury food trade (Liang et al., 2004). To this end, it is 
now listed as critically endangered both in the Chinese 
Red Book of Amphibians and Reptiles and in the IUCN 
Red List of Threatened Species (Jiang et al., 2016). Its 
current distribution is highly fragmented and several 
populations are already extirpated (Wang et al., 2004). 
Nearly all wild individuals are dwelling in mountain 
caves where they breed (Luo et al., 2009; Su et al., 2009). 
Larvae produced inside the caves and flushed out in late 
winter and early spring might mature into adults in rivers; 
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however, they have been heavily poached at the exit of 
caves with mist nets or electric fishing device in the past 
few years (Wang and Zhang, 2012). In addition, human- 
facilitated introductions but without pre-genetic screening 
also pose a genetic pollution threat (Murphy et al., 2000). 

For obligate cave dwellers, caves present formidable 
barriers to successful colonization and to movement 
between caves (Culver et al., 2006). Most cave 
populations of CGS may also be isolated due to the 
barriers or threats to immigration, including dam 
construction, alteration of river courses, pollutions, 
and poaching, and are thus predicted to have risks of 
extinction under fluctuant circumstances due to genetic 
problems such as genetic drift and inbreeding depression 
(Frankham, 2005). 

Genetic monitoring is costing less and being more 
sensitive and reliable than traditional monitoring 
approaches for conservation and management (Schwartz 
et al., 2007). A recent reduction in abundance generates 
specific patterns of nucleotide variation over time 
compared with mutation drift equilibrium. Effective 
population size (Ne), that refers to the size of a theoretical, 
ideal population with the same genetic characteristics 
as the actual population, can be assessed on the basis of 
maximum likelihood estimations or Bayesian analysis 
(Wilson et al., 2003; Kuhner, 2006). In addition, it is 
possible to detect recent bottlenecks based on detecting 
alleles’ deviations from expectations under mutation drift 
equilibrium (Peery et al., 2012; Hoban et al., 2013). 

Amphibian larvae in species where they are more 
accessible and abundant than adults are often sampled 
for population genetic study (e.g. Cosentino et al., 2012; 
Peterman et al., 2014). Although ideally CGS adults 
would have been sampled, low encounter rates of this rare 
species due to long-term poaching and lots of complicated 
shelters meant that this was impossible in karst habitats. 
Here, we identified full siblings of CGS larvae that 
flushed out of caves, addressed the effective population 
size and genetic structure of populations inside, and 
highlighted the urgency of conservation in subterranean 
habitats for this amphibian species most in need of 
conservation action, identified by the EDGE of Existence 
programme of the Zoological Society of London (Isaac 
et al., 2012). 


2. Materials and Methods 
2.1 Study Sites and sampling The three sampled 


caves are illustrated in Figure 1. Only about 300-600 
larvae were flushed out of the Shaliangzi cave (SLZ; 
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Zhouzhi County, Shaanxi province), as described by 
Wang and Zhang (2012), in January-March of 2011 and 
2012. The larvae that flushed out of the Taoyuan cave 
(TY; Nanjiang County, Sichuan province), with four 
exits at the same bank within a distance of 600 m, were 
around one thousand each year (local poachers, personal 
communication). The larvae that flushed out of the Yanxia 
cave (YX; Yanxia County, Guizhou province), located at 
the bank of a small stream, averaged 514 each year during 
2004—2006 (Su et al., 2009). 

Tail tissue of 61 and 56 injured larvae that flushed 
out of the SLZ cave in 2011 (grouped as S11) and 2012 
(grouped as S12), respectively, were provided by the 
local Fishery Bureau. Oral swabs were collected from 
48 larvae collected by villagers from the TY cave during 
2011—2012, and from 24 juveniles collected as larvae 
by farmers from the YX cave in the spring of 2008, 
following the procedure of Poschadel and Möller (2004). 
Tail tissue and oral swab samples were stored individually 
in 95% ethanol and frozen at —20°C prior to laboratory 
analysis. 


2.2 DNA isolation and characterization Total genomic 
DNA was extracted using the EasyPure'™ Genomic DNA 
kit (Beijing TransGen Biotech Co., Ltd.). The complete 
mitochondrial control region (i.e. D-loop, 770—771 bp) 
was amplified following the procedure of Tao et al. 
(2005). The PCR products were sequenced on ABI 3730 
in both directions. 

Variability of 20 microsatellites, including GS66, 
GS78, апа GS132 locus (Meng et al., 2008), CIBad01- 
CIBad16 loci (Wang et al., 2014) апа CIBad17 locus 
(with primers 5'-TTCTCACAGATGGAATCACA-3' 
апа 5'-GAAGCAGCTATGCCGTAT-3'), were detected 
following the procedure of Wang et al. (2014). 


2.3 Genetic diversity D-loop sequences were aligned 
in Mega v6.0 (Tamura et al., 2013) and calculated the 
haplotype diversity (hd) and nucleotide diversity (л 
and Ө) in DnaSP v5.10.1 (Librado and Rozas, 2009). In 
DnaSP v5.10.1, we also analyze D-loop haplotypes with 
the Tajima’s D test, which is based on the comparison 
between the average pairwise difference and the number 
of polymorphic sites. If equilibrium is not reached after a 
demographic event, negative D values are expected under 
population expansion and positive D values are expected 
under population decline (Tajima, 1989). 

Potential microsatellite genotyping errors were 
detected using Microchecker software v2.2.3 according to 
Brookfield’s approach (Van Oosterhout et al., 2004). We 
estimated allelic diversity, overall observed and expected 
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heterozygosity, and tested for linkage disequilibrium 
between all possible locus pairs using Fisher’s exact tests 
in GENEPOP 4.0 (Rousset, 2008). In GENEPOP, we also 
tested for deviations from Hardy-Weinberg frequencies, 
reflected as whether the inbreeding coefficient (f; Weir 
and Cockerham, 1984) is significantly different from 
zero; exact P values are calculated with the Markov chain 
method (dememorization 10 000, batches 100, iterations 
per batch 5000; Guo and Thomson, 1992). We performed 
a sequential Bonferroni correction to minimize type-I 
errors by adjusting P- values to correct for multiple table- 
wide tests (Rice, 1989). 

We calculated the number of rare alleles (Ar) 
(defined as alleles with a frequency <1% in all 189 
samples), and the frequencies of private alleles in each 
population. Allelic richness, which could be considered 
as an alternative criterion to measure genetic diversity 
(Rodrigafiez et al., 2008), was standardized for variation 
in sample size and calculated based on the minimal 
sample size of 24 in this study. 
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2.4 Number of breeders and effective population size 
To provide an indication of population size inside the 
karst cave, sibships of larvae and Ne were reconstructed 
using the COLONY program (Jones and Wang, 2010), 
with sets of full-pedigree likelihood method, polygamous 
mating system of both sexes, no sibship size prior, allelic 
dropout rate of 0.01 and false allele rate of 0.01 for each 
locus. The same dataset was processed 8—16 runs with 
different random seeds to confirm the reliability of results. 


2.5 Bottleneck tests Microsatellite data was analyzed 
with the heterozygosity excess test (Luikart et al., 1998) 
and M-ratio test (Garza and Williamson, 2001). We used 
the Wilcoxon two-tailed test in BOTTLENECK (Piry 
et al., 1999), as it is more powerful than the sign test 
when average sample size per site is 30 or less, and the 
standardized test requires a minimum of 20 loci (Luikar 
et al., 1998). We chose a two-phase mutation model 
(TPM), which is modified from the strict stepwise 
mutation model (SMM), because it allows for (a user- 
defined percentage of) multi-step mutations and 
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Figure 1 Position of the three sampled karst caves. SLZ: Shangliangzi, Zhouzhi County, Shaanxi Province; TY: Taoyuan, Nanjiang County, 


Sichuan Province; YX: Yanxia, Guizhou County, Guizhou Province. 
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consequently to be best fit for most microsatellite datasets 
(Piry et al., 1999). In BOTTLENECK, we chose a 5% 
multi-step mutations (Piry et al., 1999; Storfer et al., 
2014). 

The М-гайо test is based on the frequency distribution 
of allelic sizes, which is expected to have gaps after 
a bottleneck due to stochastic loss of rare alleles. The 
М-тайо is computed in each data set as simply the ratio 
of the number of alleles per locus divided by allelic size 
range for that locus. Evidence of deviation from the null 
hypothesis of demographic stability can be concluded 
if the observed value is lower than a simple threshold 
value of 0.68 (Garza and Williamson, 2001), which is 
widely used as a “rule of thumb” in conservation genetics. 


2.6 Population subdivision D-loop haplotypes in the 
three populations, together with other sequences of 
Cryptobranchidae species (GenBank code GQ368662, 
Zhang and Wake 2009; AB445787-99, AB445800-03, 
Matsui et al., 2008) were used to reconstruct the 
gene genealogy using maximum likelihood (ML) as 
implemented in PhyML 3.0 (Guindon et al., 2010). We 
use Cryptobranchus alleganiensis as an outgroup and 
conducted bootstrap analysis using 1000 ML bootstrap 
pseudoreplicates of the dataset, fixing model parameters 
at maximum likelihood values and estimating tree 
topology. 

Clusters of individuals in the sampled populations 
were also inferred using microsatellite genotype data in 
STRUCTURE 2.3.4 (Pritchard et al., 2000). Simulations 
were run using a burn-in period of 10 000 sweeps 
followed by 100 000 MCMC iterations. The number 
of possible clusters (К) was set from 1 to 12 (i.e. gene 
pools assumed), and 30 independent runs were repeated 
to guarantee convergence. Twenty independent runs with 
the highest log likelihood [Ln Pr(X/K)] of the posterior 
probability were selected to estimate the true K (Evanno 
et al., 2005). 


3. Results 


3.1 Genetic diversity Two haplotypes with one 
transversion site were found in the D-loop region in the 
SLZ population, with the haplotype diversity (hd) and 
nucleotide diversity (л, 0,,) of 0.503, 0.00053 and 0.00022 
in 2011, and 0.503, 0.00053 and 0.00023 in 2012, 
respectively. In contrast, one unique haplotype was found 
in the D-loop region in either TY or YX populations. 
Micro-Checker revealed no evidence for null alleles 
of all microsatellite loci except GS132 and GS78 due to 
stuttering. The numbers of loci that deviated from Hardy- 
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Weinberg equilibrium (P < 0.05) were 15 (S11), 16 (S12), 
11 (TY), and 4 (YX), respectively (Table 1). Mean alleles 
in the populations were 6.8 (4—8, S11), 4.9 (2—7, S12), 
4.9 (1-9, TY), and 3.0 (1-5, YX), respectively, with the 
corresponding allelic richness of 2.7 + 0.1, 2.1 + 0.1, 2.4 


+ 0.2, and 3.0 + 0.3, respectively. 


3.2 Population size Based on the best sibship 
configuration in COLONY, N, with 95% confidence 
interval was estimated to be 9 (5—24), 25 (15—44) and 17 
(8—35) in the SLZ, TY, and YX caves, respectively (Table 
2). All the 61 larvae that sampled in the SLZ cave during 
2011 were assigned into 6 families (clusters) that come 
from 7 fathers and 7 mothers. Similarly, the 56 larvae 
that sampled in this cave during 2012 were assigned into 
only 2 families that come from 9 fathers and 8 mothers. 
Combined these two years, all the 117 larvae were 
estimated to come from 11 fathers and 10 mothers. The 
breeders in the TY and YX caves were estimated to be 29 
and 15 breeders, respectively. 


3.3 Bottleneck and genetic drift The Tajima’s D values 
(P) of D-Loop haplotypes in SLZ were 2.35 and 2.28 in 
2011 and 2012, respectively. In the SLZ population, all 
microsatellite loci except CIBad06 and CIBad08 reduced 
alleles from 2011 to 2012, with the allelic richness 
reduction of 22% (2.72.1). Only two loci (CIBad05-06) 
had one new allele in 2012. Significant heterozygosity 
excess (P < 0.05) was detected at 4, 5, 9, and 2 loci in 
S11, $12, TY, and YX, with the Wilcoxon p values of 
0.648, 0.003, 0.000, and 0.196, respectively, indicating 
recent bottlenecks in these cave populations (Table 3). 
All М values at each locus in the three populations were 
lower than the threshold value (0.054—0.429 vs. 0.68, 
Table 3), also indicating serious population decline 
recently. 


3.4 Genetic subdivision There were 16 variable 
nucleotide sites across the D-loop region for the 
three populations and the overall p-distance between 
populations is 0.011. The rooted ML phylogeny reveals 
that mtDNA haplotypes of CGS are separated into three 
reciprocally monophyletic clades (Figure 2). The TY 
population was closer to the SLZ population but further 
to the YX population (Figure 2). Compared with the 
D-loop haplotype in YX population, all haplotypes in 
SLZ and TY populations have one deletion A (i.e. 770 
bp). STRUCTURE analysis revealed one peak at K=3 of 
AK value for the three Andrias populations (Figure 3a). 
When three clusters were assigned, all individuals have 
the admixture proportion of > 82% to be a cluster and no 
mixture between different populations (Figure 3b). 
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Figure 2 Midpoint-rooted maximum likelihood (ML) phylogeny of Chinese giant salamander mtDNA D-loop haplotypes. The numbers on 
branches represent ML bootstrap values (1000 replicates). Cryptobranchus alleganiensis are used as an outgroup. 


4. Discussions 


4.1 Genetic variance and population subdivision 
Only one or two D-loop haplotypes were detected in 
the three sampled populations. Similarly, only one to 
three haplotypes of Cytb or D-loop were found in each 
population of Japanese giant salamander (Matsui et 
al., 2008) and North American hellbender (Sabatino 
and Routman, 2009). However, at least 27 haplotypes 
of D-loop were detected in CGS until now (Tao et al. 
2005), indicating that the mitochondrial genetic variation 
of Chinese giant salamander mainly happens between 
geographic units but not within populations. 

However, a much higher haplotype diversity (hd, 
0.842—1.000 vs. 0.503) and nucleotide diversity (0.00496— 
0.01114 vs. 0.00065) were found in previous studies (Tao 
et al. 2005; Fang et al. 2008; Table 4), as their samples 
were from a much wider region (such as the upper 
Yangtze River; Tao et al., 2005) or from farms, where 
adults from far different rivers were stocked together. 

Previous studies showed genetic divergence of CGS 
in China even though no distinct geographic structure 
were found partially due to human relocation (Murphy 
et al., 2000; Tao et al., 2005; Tao et al., 2006; Yang et al., 


2011). Only one D-Loop haplotype was found in the TY 
and YX population, and two haplotypes differentiated 
with only one nucleotide site were found in the SLZ 
population, which may indicate that no introduced 
genome existed in our cave populations at present. 
However, the clear genetic divergence among the three 
cave populations in terms of mitochondrial (Figure 2) 
and nuclear polymorphism (Figure 3) further support the 
downstream associated detrimental genetic pollution due 
to governmental-oriented releasing activates. 


4.2 Small population size and genetic drift Based 
on some sporadically published breeding parameters, 
including brood amount of 400 eggs/kg for mother, egg 
fertility rate of 70%, and embryonic hatchability of 5%, 
Su et al. (2009) estimated that the sum weight of female 
breeders were 37 kg inside the YX cave. Considering 
that the female breeders (0.35—3.5 kg; Wang et al., 2000) 
may reach large size (e.g. 1.5 kg) in caves because of the 
absence of poaching, the estimation of Su et al. (2009) 
were three times to our sibship analysis (24 vs. 7 mothers, 
Table 2). One reason may be that their hatchability was 
underestimated (e.g. > 10% of A. japonicus, Suzuki et al., 
2000). 

Populations with a large Ne are resistant to loss of 
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Table 2 Sibship clusters and effective population size (Ne) analysis using COLONY program based on the microsatellite genotypes of larvae 


from SLZ, TY and YX. 


SLZ 
TY YX 
2011 2012 2011+2012 

Sample size 61 56 117 48 24 
No. of fathers К 9 11 17 8 
No. of mothers 7 8 10 12 7 
Total parents 14 17 21 29 15 
Clusters 6 2 6 1 1 
Ne (95% СТ) 9 (5—24) 8 (4-23) 9 (5—24) 25 (15—44) 17 (8—35) 


Table 3 Results from BOTTLENECK (ТРМ model) and estimates of М (the ratio of number of alleles per locus divided by allelic size rage 
for that locus), as an estimate of putative bottlenecks for populations of SLZ (S11, S12: п = 61 and 56 respectively in 2011 and 2012), TY (л 


= 48), and YX (л = 24). 


cae Probability under TPM M 

511 512 ТҮ. ҮХ 511 512 TY YX 
9566 0.339 0.046 0.013 0.263 0.156 0.133 0.333 0.066 
9578 0.419 0.043 0.365 0.330 0.276 0.280 0.308 0.333 
98132 0.016 0.454 0.062 0.264 0.226 0.161 0.429 0.308 
CIBad01 0.410 0.319 0.004 0.006 0.207 0.172 0.280 0.122 
CIBad02 0.332 0.233 0.056 0.236 0.212 0.152 0.235 0.154 
CIBad03 0.201 0.220 0.123 0.425 0.241 0.235 0.171 0.190 
CIBad04 0.204 0.337 0.389 0.462 0.156 0.235 0.308 0.069 
CIBad05 0.162 0.212 0.213 0.320 0.242 0.135 0.333 0.103 
CIBad06 0.012 0.034 0.001 0.463 0.212 0.189 0.220 0.190 
CIBad07 0.354 0.072 0.068 0.103 0.171 0.190 0.143 0.154 
CIBad08 0.028 0.298 0.005 0.294 0.200 0.162 -- 
CIBad09 0.149 0.013 0.137 0.038 0.171 0.152 0.200 0.222 
CIBad10 0.145 0.001 0.116 0.178 0.133 -- 0.118 
CIBad11 0.013 0.301 0.033 0.162 0.240 0.222 0.171 0.098 
CIBad12 0.360 0.174 0.171 0.191 0.212 0.172 0.138 0.176 
CIBad13 0.427 0.074 0.020 0.190 0.231 0.333 -- 
CIBad14 0.143 0.166 0.021 0.145 0.163 0.102 0.240 0.190 
CIBad15 0.057 0.128 0.022 0.103 0.113 0.152 0.333 0.222 
CIBad16 0.124 0.279 0.359 0.068 0.182 0.190 0.146 0.111 
CIBad17 0.221 0.269 0.028 0.292 0.242 0.207 0.170 0.054 
Wilcoxon (P) 0.648 0.003 0.000 0.196 


genetic diversity through genetic drift and inbreeding. 
It is estimated that the minimum viable Ne should not 
be less than 50 in the short term and should be greater 
than 500 in the long term to maintain genetic diversity 
(Jamieson and Allendorf, 2012). If we use this as a guide, 
it is likely that the sampled cave populations (with Ne 
< 44) will decline in genetic diversity rapidly. In the 


SLZ population, we found that Ne reduced from 9 to 
8 and the allelic richness declined at a faster rate than 
heterozygosity (2.7—2.1 vs. 0.710—0.671), primarily 
due to the loss of rare alleles (2.25—0.60; Table 1). 
More alleles (1—5) deviated from the mutation-drift 
equilibrium also showed that this population decreased 
to a threshold and face genetic bottlenecks. Arguably it is 
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Figure 3 Clustering assignment depending on the Bayesian method under an admixture model obtained by STRUCTURE software. (a) 
AK values as a function of K based on 20 runs, indicating the most likely number of three genetic clusters, (b) Output of genetic clusters 
identified, each individual is represented by a single column that is divided into segments whose size and color correspond to the relative 
proportion of the animal genome corresponding to a particular cluster. Populations are separated by black lines.. 


allelic richness rather than heterozygosity that reflects the 
long-term evolutionary potential of a population, since 
it is the number of genetic variants in a population that 
determines the material available for natural selection to 
act upon (Tracy et al., 2011; Weiser et al., 2013). Recent 
meta-analyses have shown that the loss of diversity at 
neutral loci following a population bottleneck is often 
accompanied by an equivalent or greater loss of diversity 
at other highly variable adaptive genetic loci, such as 
major histocompatibility complex (MHC) genes that are 
important in disease resistance (Radwan et al., 2010; 
Sutton et al., 2011). Field surveys indicated that the 
larvae production of wild populations decreased sharply 
throughout its range and even disappeared at some sites 
in recent years (Luo et al., 2009; Su et al., 2009; Wang 
and Zhang, 2012). If this is the case, our results suggest 
managers to make concrete efforts to maintain allelic 
richness and enlarge effective population size for the 
long-term viability of cave populations. 


4.3 Conservation suggestions Some reservoirs or rivers 
were reported to find large adults in recent years, but such 
places, polluted with domestic sewage or agricultural 
chemicals or out of this species’ range (such as Beijing 
and Shandong Province), were usually unfit for breeding 
and seemed to harbor released or escaped captive-bred 
individuals. The natural breeding caves like SYZ, TY, and 
YX provide very important refuges to CGS. The larvae 
that flushed out caves are important recruitment to cave 


populations and need to keep away from poaching for 
farming stock. 

Subterranean ecosystems generally exhibit little or no 
primary productivity and rely almost entirely on resources 
produced outside caves (Culver and Pipan, 2009). No 
larvae were flushed out of the YX cave after the water 
inlet was blocked with concretes to prevent poisoning in 
2008, which obviously changed the inner hydrological 
conditions and cut off the food resources of breeders 
inside. Floods are important sources of organic input into 
cave systems, however, it may reduce ovipositon sites 
or cause asphyxiation of embryos. These could be the 
possible hypotheses why larvae were not observed out of 
the SLZ cave in 2013. 

Land use decisions, such as unsanctioned quarrying, 
alteration of hydrology, illegal logging, and increased 
tourism, can have immediate and serious impacts on karst 
terrains (van Beynen et al., 2012). A number of limestone 
cave-dwelling amphibians, such as Red-spotted toothed 
toad (Oreolalax rhodostigmatus) and Wuchuan odorous 
frog (Odorrana wuchuanensis) listed as vulnerable or 
critically endangered species by IUCN, are particularly 
sensible to such human activities. 

To aid the recovery of wild population, the Chinese 
Ministry of Agriculture organized a series of releasing 
activates of CGS from 2002 and about ninety-two 
thousand artificially bred juveniles were introduced into 
the wild during 2010-2013. However, most releases 
of CGS failed to boost the wild populations due to the 
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limitation of factors such as water quality, food supply, 
and safety (Luo ег al., 2009). We suggest that sustainable 
forestry, landfills preventing groundwater contamination, 
and nature reserves extension to connect breeding caves 
instead of releasing may be the sole measure to save this 
flagship species in karst habitats. 
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